

effects = c(-4, -3, -2, -1, 0, 1, 2, 3, 4)

mother = c(0, 0, 0, 1, 1, 1, 2, 2, 2)
father = c(0, 1, 2, 0, 1, 2, 0, 1, 2)


child = c(0, 0, 1, 0, 1, 2, 1, 1, 2)

childBV = sum(effects*child)
childPheno = childBV + rnorm(1, 0, 1)

evaluate = function(genotype, childPheno){
    estimate = sum(genotype*effects)

    prob = dnorm(childPheno, estimate, 1)
    return(prob)
}

getPossibleGenotypes = function(g1, g2){
    if(g1 == 0 & g2 == 0) return(c(0))
    if(g1 == 1 & g2 == 0) return(c(0, 1))
    if(g1 == 2 & g2 == 0) return(c(1))
    if(g1 == 0 & g2 == 1) return(c(0, 1))
    if(g1 == 1 & g2 == 1) return(c(0, 1, 1, 2)) #To account for the two ways of getting a 1 state.
    if(g1 == 2 & g2 == 1) return(c(1, 2))
    if(g1 == 0 & g2 == 2) return(c(1))
    if(g1 == 1 & g2 == 2) return(c(1, 2))
    if(g1 == 2 & g2 == 2) return(c(2))

}

getPenetrance = function(index, options, probs){
    geno = c(0, 1, 1, 2)
    pen = c(0, 0, 0, 0)
    for(i in 1:4){
        pen[i] = sum(probs[options[,index] == geno[i]])
    }
    pen = pen/sum(pen)
    return(pen)
}
evaluatePermutations = function(mother, father){

    vals = list()
    for(i in 1:length(mother)){
        vals[[i]] = getPossibleGenotypes(mother[i], father[i])
    }

    options = expand.grid(vals)

    probs = apply(options, FUN = evaluate, MARGIN = 1, childPheno = childPheno)

    penetrance = lapply(1:length(mother), getPenetrance, options, probs)
    return(do.call("cbind", penetrance))
}



